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Section 1 


INTRODUCTION 


This is the second of a series of reports on the develop- 
ment of finite-element analysis procedures for computation of 
linear elastic stress intensity factors associated with cracks 
in common aircraft structural details. Finite-element analysis 
serves the important purpose of providing Ky and Kir solutions 
for situations in which she boundary geometries are too complicated 
to permit convenient treatment by the classical methods. 
The first report in the series | 1] summarized the formula- 
tion of two basic building blocks which are used in the analyses: 
the well-known bilinear isoparametric quadrilateral and the PCRK59 
-assumed-stress hybrid crack-containing elements. The former element 
- has been widely accepted and appears in most general-purpose finite- 
element programs. The latter element was originally developed at 
ASRI. in 1972 [2], and has subsequently been put through extensive 
‘per cosnanke tests to characterize its behavior as a function of 
_ shape distort ons [ 3}. 
Earlier work (4, 5] has demonstrated that the assumed- 
stress hybrid method permits accurate computation of stress 
intensity factors, while preserving Matrix Displacement Method 
conventions in the glohal programming environment (data input, 


assenbly, matrix solution erocedure, etc.). In the present work, 


the hybrid method 1s used auain to formulate another special- 


purpose element for modelling the region around a fastener hole. 
The new element is combined with the quadrilateral and PCRK59 to 
represent a skin tension panel containing an open fastener hole. 
The ASRL FEABL program, 6) is again used for global analysis. 


Additional features belonging to FEABL which have not yet been 


outlined briefly in this report. These features pertain to the 


technique of substructurir bh, a Gauss elimination algorithm. 
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formally documented are employed in the preseit analysis, and are 
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Section 2 


NEW DEVELOPMENTS 


Modularization of the finite-element model of a structure 
becomes as important as modularization of the basic program when 
the application is computation cf stress intensity factors. The 
objective is to reduce unnecessary arithmetic cperations and 
computer printout to a minimum. In particular, the stresses and 
displacements are of little interest (and of no interest at all 
at any distance away from the crack tip) except for verifications 
of accuracy during development. The substructuring technique, 
used for many years in airframe analysis programs*, provides the 
required modularity. Subroutines for substructuring were developed 
as add-ons to FEABL as pirt of an unrelated project** and have 
been used in the present work. The substructuring algorithm ard 
the functions and conventions of the add-on subroutines are dis- 
cussed briefly in this section to provide a complete picture of 


the panel analysis program. 


The so-called “near-field” region around a fastener hole 


forms one such substructure module. The first attempt to model 


——_ 


*For example, the Boeing ASTRA program. 


**Sponsored by the Xerox Corporation (Dr. T. C. Soong, technical 


this region with conventional quadrilateral elements resulted in 
| 
monitor) . \ 

t 
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severe mesh distortion, numerous unwanted degrees of freedom, 


and solution accuracy problems. In consequence, the hybrid 
method was called upon, resulting in the HOLEI special-purpose 
element for this module. The history of these investigations is 


traced in Subsections 2.3 and 2.4. 


2.1 Substructuring 

Figure 1 illustrates the finite-element model of a square 
panel which contains a single fastener hole offset from the 
panel centerJine. A detail in the lower part of the figure shows 
how a PCRK53 element might be inserted to represent a small crack 
emanating from the hole. It is evident that the complete global 
solution (displacements at each of the many nodes in the model) 
constitutes unwanted informaticn if the analysis ic to be repeated 
while the crack location is varied parametrically around the 
fastener hole. 

A computationally more efficient procedure is obtained by 
subdividing the model into far-7ield and near-field substructures, 
as shown in Fig. 2. ‘The broken line in Fig. 2 corresponds to 
the edges ABCD in Fig. 1. If the element stiffnesses, prescribed 
nodal forces and nodal displacements are first assembled for the 
far-field substructure, the resulting equation system may be 
partially solved. In effect, all of the information pertaining 
to the far field substructure is transterred to th: inter-sub- 
Structure boundary ABCD. This information may then be treated 


as a “super-element", which can be assembled with a succession 


of near-field substructures containing a crack in different loca~ 
tions. As a result, many fewer degrees of freedom are processed 
in the parametric stress intensity computations. 

The substructuring process may be represented formally as 
follows. Let the assembled equation system Kq = Q for a sub- 


structure be partitioned into: 
K., Kis 


(1) 
Ko Kus 4s 


> > 


stiffness coefficients 


in 
iT] 


it 


(unknown) nodal displacements 


prescribed nodal forces 


1O>4aQ 


ca] 


and where K,, = K Subscripts I and B refer to nodes and degrees 


BI ~IB° 
of freedom which are considered respectively as “interior" (to 
be eliminated) and “boundary” (to be retained for subsequent 
assembly along the inter-substructure boundary). The interior 
degrees of freedom are formally eliminated by solving the first 


of Eqs. 1 and substituting into the second; 


wd A on" ) A 
a8. Keb) Bute Be 


ATER STE oe THES 


the unknowns Gy pe In practice, the process may he programmed 


to eliminate one degree of freedom (one equation) from the system 
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at a time, so that inversion of Kir is avoided. The entire 


algorithm then consists of as many passes through the assembled 


equations as there are degrees of freedom to be eliminated, a 


oy 
| 
ee: 
a A rearrangement of the terms in Eq. 3 leads immediately to: | 
4 K? Q° ‘ : 
E ~st,” “ap (4) i 
4 where ae of) are the statically condensed stiffnesses and 
3 forces, given by: 
3 (c) 7 
q Kee" Kae K Kas Big (5) | 
: Ate) A ~iLia ‘ 
: Q,= Q,- K,,;,@ ) 
By ~ w~Biwlt OF : 
q ae : (6) : 
4 It is apparent from Eqs. 5 and 6 that the substructuring ; 
‘ process may be realized as a computing algorithm independent of i 
B x 
‘ ; 


proceuure quite similar to the Gauss-Cholesky factoring methods 


normally used to solve the whole equation system [6 ]. 
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e Substructuring places an additional burden on the program 
4 user when the global software stores K as a band~matrix. It can 


be shown by tracing the details that the elimination algorithm 


inflates the stiffness matrix bandwidth if dy and dp are arbitrarily 


) 
Ri 

&) 

“A 

Z 

S 

5 

.¢ 
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interspersed. However, if qy always appear as the first degrees 
of freedom in the global numbering sequence, the bandwidth of K 
is not affected by elimination. This restriction has been placed 
on the FEABL add-on subroutines, and it requires some care in the 
choice of global numbering sequence and substructure boundary 


locations on the user's part. 


The foregoing derivation of the substructuring algorithm 
assumes that all the prescribed quantities are nodal forces Q, 
while all of the displacements q are unknowns. However, if 
restraints are appiied to the substructure such that a subset 
qr of the interior degrees of freedom qd; are prescribed, Eqs. 2 
ehrougii 6 may still be used. The only difference is that the 
rows and columns of K corresponding to ay are decoupled from 
the equation system*, while the right hand side of Eqs. | is 


replaced by: 


Q,} =7Q oa Q- K 4s, 5 (7) 


° a” 


2.2 Implementation of Substructuring in FEABL 


Three add-on subroutines have been programmed and verified 
for the gi} -al computation tasks associated with substructuring. 
Subroutine STACON executes the Gauss elimination algorithm on a 


band-matrix which has been assembled and stored in accordance with 


*Sce Ref. 6, Subsection 3.2.5. 


the existing FEABL conventions. A normal MAIN program is written 
to control the FEABL software up to the point at which K would 
be factored in a conventional analysis. However, instead of 
factoring the instructions: 

ISGN = 

CALL STACON (ISGN,NID,RNAME, INAME) 
is programmed, where 


ISGN 


i] 


Matrix condition contro. parameter 
NID = Total number of interior degrees of freedom 
RNAME, INAME = Real (floating point) and integer 
(€ixed point) variable names assigned 
to the FEABL DATA vector 
The control parameter ISGN must be set to +1, -] or 0 before STACON 
is called. A value of +1 aborts the run if K is not a positive- 
definite matrix. A value of -1 aborts the run if K is singular. 
A value of 0 always results in return to the MAIN program with 
ISGN reset according to the actual matrix condition encountered: 
positive-definite (+1), non-positive but nonsingular (negative 
integer), or singular (0). The positive-definite cption is used 
in conventional finite-element analyses. The nonsingular option 
is used in special cases, e.g., Hermann's principle for which 
ideally incompressible elements include volumetric constraints 
in the assembled equation system. The option ISGN=0 may be used 
generally, in programs for which the user wishes to retain control 
for debugging or other purposes when an error condition is encountered. | 


Subroutine STACON carries out the elimination process in place, 
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as shown schematically in Fig. 3. After execution, the statically 
condensed force vector appears in INAME(IQ+NID) to INAME(LQ) aad 
the statically condensed stiffness coefficients appear in INAME (J) 


to INAME(LK), where: 


J=IK+INAME (IKGUNT+NID) +NID+1 (8) * 
Tne shaded areas shown in Fig. 3 contain the back-substitution 
information Ke Q, and ae Kips which may be used to solve for 
Qy (Eq. 2} after qn has been computed. 

The statically condensed stiffnesses and forces may simply 
be transferred to another storage area and treated thereafter as 
a “super-element" if desired. A qeneral programming sequence to 
execute the transfer is given in Appendix A. If this is done, 
the substructure may later be assembled into a complete structure 
with subroutine ASEMBL or subroutine ASMLTV. A list of local-to- 
global degree-of-freedom connections must be given for the “super~- 
element", just as is done for any ordinary element. However, 
while the numbering convention for an ordinary element is fixed 
by the manner in which the element subroutine has been programmed, 
the convention for a “super-eiement" is derermined by the user's 


sequence of global numbering for assembly of the substructure. 


For example, if the far-field substructure in Fig. 2 has been 


*IKZUNT is another address control parameter. See Ref. 6, 
subsections 2.1.3, 2.2.2 and 2.2.6. 
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numbered such that the boundary degrees of freedom are globally 
assigned in the order A+B-C+D+A, then the connections of the far- 
field “super-element" to the near-field portion of the model must 
be given in that order. 

A second option available to the user is to leave the statically 
condensed substructure in place and to employ a second DATA vector 
for assembly and solution of the complete structure. Add-on sub- 
routine ASMSUB has been programmed for this option. Subroutine 
ASMSUB performs a function similar to the conventional assembly 
routines (ASEMBL,ASMLTV), but is capable of assembling ane on 
directly from one DATA vector into another. Since two DATA vectors 


must be manipulated simultaneously, extra address control parameters 


are required. The beginning of the program might appear as follows: 


DIMENSI@N RNAME] (xxx), INAME] (xxx), RNAME2(yyy), INAME2 (yyy) 
EQUIVALENCE (RNAME1(1), INAME1(1)) 

EQUIVALENCE (RNAME2(1), INAME2(1)) 

COMMBN /SIZE/ NET, NDT 

COMMON /BEGIN/ ‘CON, IKOUNT,ILNZ, IMASTR, IQ, IK 

COMNGN /END/ LCJN, LKSUNT, LENZ, LMASTR,LOQ, LK 

COMMZN /SIZESS/ NETSS,NOTSS ,NIDSS 

COMMON /BEGSS/ IBEG(6) 

COMMZN /ENDSS/ LEND (6) 


. 
» 


s 


Now suppose that a substructure is assembied and statically 
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Se 


| 
} 
1 
| 


condensed in vector (RNAME1, INAME]): 


e 
e 


ISGN=1 

CALL STACON (ISGN,NID, RNAME], INAME1) 
At this point, the substructure's control information must be 
cleared from the labelled COMMON areas to permit formation of 
the complete structure. The information is transferred to the 
extra areas: 

NETSS=NET 

NDTSS=NDT 

NIDSS=NID 

IBEG(1)=ICON 


LEND (6) =LK 


° 
. 


The program continues with a conventional generation of size, 
connections, assembly, etc. for the complete structure in vector 
(RNAME2, INAME2) until the substructure is to be assembled. At 
this point, the user simply programs: 


CALL ASMSUB(LNUM, RNANE1, INAME], RNAME2, INANE2) 


° 
° 


* 


where LNUM represents tho element number assigned to the sub- 
structure, considered as a “super-element" in the complete 
structure. 

There are some cases for which the solution Q, may be 
required in a particular substructure, Subroutine QBACK has 
been proqrammed as the third add-on to execute the required back- 
substitution solution for qy° Subroutine QBACK is designed to 
work with ASMSUB in the multiple DATA vector environment. Suppose 
that the above example is continued to the point ut which assembly 
and constraint of the structure equation system are complete. 


Factoring and global solution now follow: 


. 
. 


. 


CALL FACT(ISGN, RNAME2, INAME2) 

CALL SIMULO(ENERGY, RNAME2, INAME?) 
At this point, vecter (RNAME2, INAME2) contains the solution gq 
for the complete structure. Part of q is, in fact, the boundary 
displacement solution Yn for the substructure. An interior solu- 
tion may be obtained by continuing with: 


CALL QBACK (LNUM, RNAME], INARE1, RNAME2?, INAME2) 


After execution of QBACK, the substructure qlobal solution iq, ig.) 


appears in vector (RNAME],INAME1]). The solution will also be 


le 


printed by subroutine QBACK. 


2.3 Experiment with Conventional Mesh 


The first attempt to model the panel and fastener hole was 


made with conventional quadrilaterals and a S-node hybrid mesh- 
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} 
| 
a expander element. Experience with analysis of an attachment lug 
detail {1 } indicated that at least 24 guadrilaterals shoul’ be 
Placed around the hole to obtair accurate solutions. One then 
faces the following choices for the panel: 
| 1. Mesh expansion in the near-field region to reduce the 
| number of divisions below 24 on the near/far-field 
| boundary (ABCD in Figs. 1 and 2). 
2. Acceptance of 24 divisions on the near/far-field Soundary 
and continuation of tais detail into the far-field region. 
3. Acceptance of 24 divisions en the boundary, with some 
fresh expansion immediately outside the near-field reaicn, 
The first choice was judged te be unlikely to give accurste computed 
displacemants rear the hole because of the severe shape distortions 
which the mesh-expander clements would have caused. The second 
choice was deemed to be unacceptable in that too many unwanted 
degrees of freedom would be placed in the far-field region, making 
either cne-staqe solution or substructuring computationally ineffi- 
cient. The third choice was consequently selected as the rest 
compromises between computational efficiency and solution accuracy 
near the hole. 
Pagure 4 illustrates a typical mesh for a panel with the 


fastener hole centered. Extra detail is kept in the far-field 
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region to the right of the hole in order to permit the hole to 
be offset continuously toward the cight edge of the panel while 
avoiding extreme aspect ratios in the far-field elements. Tests 
with this model quickly showed that the asymmetric mesh grading 
led to highly asymmetric behavior of the computed solutions. 
Figure 5 illustrates a typical example. Computed horizontal 
displacement of the left and right edges of the panel are shown 
for the case of uniform vertical tension and the fastener hole 
centered. The inaccuracy of the model was judged to have been 
caused by the combination of asymme ric mesh grading, aspect 
ratio effects, and shape distortions in the near-field region. 
The conclusion was that conventional finite-element models 
could be made to wor for the panel structure only with very 
fine detail, and therefor at an unacceptably high computing cost. 
Lest it be thought that too much accuracy has been demanded from 
the analysis, we remind the reader that more than an engineering 
solution is required to achieve the ultimate goal. In considering 


as 


the accuracy of the K, and Ky solutions, it has been shown { 1, 31] 


TI T 
that shape distertions of the PCRKSI clement alone can account 

for as much as 3 to 5 percent error. Furthermore, direct verifica- 
tion of accuracy by csomparison with independent methods is possible 


only for K, with cracks emanating horizentally from the fastener 


i 
hole. Hence, a finite-element model capable of computing highly : 
accurate displacements near the fastener hole is required to 


achieve engineering accuracy in the final answer. As a result, 


the assumed-stress hybrid method was investigated as a possible 
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alternate way of handling the transition from the circular hole 


geometry to the cartesian geometry of the near/far-field boundary. 


2.4 Creation of a Hybrid Fastener Hole Element 


The energy principles on which hybrid elements are based 
have been reviewed previously [11. The energy principle ll, was 
chosen for the present investiyation, giving the general form of 


an element stiffness matrix as: 


DKA A ARBRE OREN LEANNA ESTAS SVP OTE LS TRUSTS MN NST be AF IE OI, 


wont F 

k=GHG : 
: ~ (9) q 
Be where ‘ 
ay ; 
: r 
: G= tt (yp) L da ao | 
Pan ~ 
: oA 
e | 
| 
H=-t{PSP4A 
E ~ A nw (11) : 
| 
re and where 
f A = Element area { 
: 
s dA = Element boundary 
f 
Ee L = Matrix of assumed-displacement interpolation 
: ra 
functicns defined on 9A. } 
& 
i N = Matrix of direction cosines for the outward 


nermal to dA. 


ne 


| 
P = Matrix of assumed-stress shape functions. 
S = Matrix of.elastic compliance constants. 
: | 
t = Elemeat thickness (assumed to be unity without 


luss of generality). 
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{wo fundamental properties af hybrid elements make them useful 
in linear elastic plane stress and plane strain analyses: 

J. Great flexibility in chvice of element shape is permitted, 
since L need be Gefined only in terms of an arc length 
coordinate from node to node along 3A, and since the 
shape functions P need only satisfy the continuum 
equilibrium equations. 

2. Many shape functions P can be chosen to satisfy equilibrium 
and at the same time allow the element to closely mimic 
a particular local stress distribution. 

Both properties have been used to advantage in the development 


of an element for the near-field region. 


Practical application of the hybrid method requires that the 


analyst account for known local »ehavior in the formulation of 


his element. In the present case, one reasonably expects to see 
sin(26) and cos(268) components in the stress distribution near 

a fastener hole in 4 panel under tension, as well as terms indepen- 
dent of @. Second, one should expect to see the near/far-field 
boundary deform in the manner shown schematically in Fig. 6. 
Simulation of this behavior requires that the element have at 

least three nodes (corners and mid-edge) along its near/far~field 
boundary edge. Third, the element must be kinematically stable, 
since it will be placed in the mesh so0 as to separate two conven~ 


tionally modelled regions. Kinematic stability is assured if: 


m2zn-3 (12) 


fer plane stress and plane strain elements, where 


m = Total number of assumed~-stress shape functions. 


n = Total number of nodal degrees of freedom (dis- 
placements) in the element. 

Finally, at least 20 nodes should be provided along the circular 
inner boundary, which may serve either as the fastener hole boundary 
or as an interface to which quadrilaterais and PCRK59 elements 

may be coupled. 

The kinematic stability condition results in a quick determina~ 
tion that modelling of the entirr: near-field region with a single 
element would be impractical. For example, assumption of the 
minimum number of nodes (8 along the outer, 20 along the inner 
boundary) results in a 56-degree of freedom element which requires 
m > 53 assumed~stress shape functions. It is seen from Eqs. 93 and 
1l that H is an mxm matrix which must be inverted to compute the 
element stiffness matrix. Furthermore, H is fully populated. The 
single-element approach was rejected in view of the difficulty 
in finding 5. shape functions and the computational inefficiency 
associated with inversion of a 53 x 53 matrix. 

A 90-degres segment of the near-field region was the first , 
geometry for which an element developmant was attempted. The 
element shape, shown in Fig. 7, encompasses one quarter of the 
circular boundary and one edge of the near/far-field boundary. 
With 20 divisions allowed around the hole, the element possesses 


6 inner nodes in addition to its 3 near/far-fileld boundary nodes. 
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The element can be subjected to rotation transformations and 
repeated assembly to model the entire near-field region. The 
kinematic stability requirement now leads to m > 18 - 3 =15 for 
the minimum number of stress shape functions. However, since 
rotation transformations are contemplated, insensitivity to 
orientation enters as another practical requirement. This leads 
to a choice of m=16, as explained in Subsection 2.5. Tests of 
the 9f-degree element highlighted a problem of inaccuracy in the 
integration for the H matrix (Eq, 11), which was judged to have 
been caused by the presence of terms in sin(46) and cos(46) in the 
stress shape functions Pp, 

The approach finally adopted divides the near-field region 
into eight 45-degree elements. Figure 8 illustrates the basic 
element, which possesses 4 nodes along its circular boundary. 

The element possesses 12 degrees of freedom, and therefore m>9 is 
required. To preserve orientation insensitivity, m=10 is chosen, 

so that all sin(26) and cos(26) terms are included in P. The entire 
near-field region is assembled by first computing ky for the basic 
element and then subjecting k) to a series of transformations for 

the seven other elements shown in Fig. 9. 

2.5 Final Formulavion of Hybrid Hole Element 

The choice of assumed-displacement interpolation functions 
for the hybrid element is dictated by .he need to maintain inter- 
element compatibility along its edges. Since the hole element 


must be capable of coupling to a quadrilaterai or a PCRK5$ segment 


between gach pair of nodes, linear interpolation must be used. 
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In general, 


fain} : [22 0 4 0 4 0 sae 41 
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tie 


where Gyr Agr +++G)o are defined in Fig. 8, and where 


u(&) = Horizontal edge displacement. 
v(4) = Vertical edge displacement. 
& = Length of edge between a pair of nodes. 


The are coordinate 4 in Eq. 13 is a relative measure of distance 
along any edge, with the positive sense taken as counterclockwise 
for integration along 3A. Note also that the nonzero terms in 


L appear in different positions for different edges. For example, 


Lb 2f[9°0% 144 0 HO ae 

~ 34 0000 0 14-44 6% go0o00 (14) 
for the edge between nodes 3 and 4, while 

E (“2 000000001-m 0 

~6i low 00000000 0 14/2 (15) 


for the edge between nodes 6 and 1. In a similar manner, 


A* 0,434 at nodes 3,4 
(16) 
4° 0,hj, at nodes 61 


and so forth. These changes of definition do not cause difficulties 
in the boundary integration (Eq. 10) because the computation is 
zarried out numerically, one edge at a time. 

Assumed-stress shape functions are selected most conveniently 
by direct adaptation from classical elasticity solutions in polar 


coordinates. Let an Airy stress function P (x, 8) be defined by: 


198 1 a ave 19F 1 9S 
tee Oe tri aet eg" are We T2360 TF Brae (17) 


Then the equations of elasticity for homogeneous isotropic material 


reduce to: 


V*V' dre) =0 (18) 


where 
22 ,lod yt 
Vie art tr ar” Ft age (19) 
By assuming as an expansion for the stress function: 
fo) otal ic ] 
Slne)= 9 Ir) + 2 f (r)L sin(ne) or cos (ne) (20) 
ne 


it can be shown that Eq. 18 reduces to a set of equidimensional 


() (yg) 


ordinary differential equations in $¢ (CL) ssc. y for 


which the following solutions are obtained [7 F: 


pr) = a,+b, mr egris dr bar (21) 


in) cn -n 2#n 2-n 
Priri= a thal tr tdyr 35 n2Zz (22) 


Stress distributions corresponding to Eqs. 21 and 22 are obtained 


by substituting back into Eqs. 17[7]: 


ft 


ae be 4 Zen t d, (142 Mar) 


= 
i] b, 2 
pe =r 2¢, +d, (3+2Lur) (22) 
vt! =O 


and -2 -N-2 
(ae = [ ay (n-n)r” 7 bn (nen?) r +, (gen-p*)r" + 


+d, (2-n-n?)r" Jf sin(né) or cosine) | 
: 22, 
ToL ay (nten)r” 24, (n*en)r +c, (aeanen*) r+ 
(24) 
td,(z-3nent) I{sentne) pr costo) | 


ta} n-z -n-2 
Cre = La, (n?-»)r -b, (wren)r tly (men) er” + 


- dy, lwt-n) ne Th sesstng) or sintno) | 


*The solution fer a) (r) is of no interest in the present work. 


Terms in sin@ and cos9 correspond to problems in which a dis-~ 
continuity has been introduced, e.g., re-welding of an annulus 
after a segment has been cut out. 
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The undetermined coefficients bo. Coy arse: qd. in Egs. 23 and 


re) 
24 play the role of assumed-stress degrees of freedom, which are 
eliminated when the element stiffness matrix is formed {1}. The 
remaining r6O-dependent portions of the terms are placed in the 
shape function matrix P. 

Shape functions have been chosen from Eqs. 23 and 24 as 
follows. For the §-independent terms, logarithmic behavior is 
not particularly desirable; only the by and Cy terms have been 
retained. This choice is admittedly arbitrary. To complete P, 
enough 6-dependent terms are chosen from 26, 40, 80,... behavior 
to satisfy kinematic stability and orientation insensitivity. 
The latter requirement demands, e.g., that all 29-dependent terms 
be retained if any are required for stability. It is apparent from 
Eqs. 24 that there are eight n6-dependent terms for each n when 
sin(n6) and cos(n@) terms are included. Thus, selection of the 
26-dependent terms provides a total of 10 shape functions, just 
enough to satisfy the requirements for the 45-degree element 
(The 90-degree element requires addition of the 46-dependent 
terms as well, making a total of 18 shape functions.). The shape 
function matrix, as finally adopted for the 45-degree element, 


is given by: 


-2:032@ -bcos20/r4 0 "4$cos20/r? -236n28 
P (68)*1 20328 = S cos2e/t* i2r*cos 20 0 2sin20 
Z23in28 -bsiniesrt brisen28 ~2sintIO/r® -200320 : 
~bsinzb/r# =O ~dsintolr? = Sfre 2 (25) ! 
bsinzo/r? tar*sinz6 0 ~UY/pd 2 
bcos26fr4 -br*1es26 — 2e0828/r* o Oo 
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The numerical integration required for computation of the 
4 matrix (Eq. 11) is accomplished most conveniently in polar 
coordinates. However, computation of G (Eq.10) actually requires 
the cartesian surface tractions which correspond to the assumed- 


displacement field: 
T - 17, TL*N fo, Syy wy} NP (26) 


where 8 represents the vector of stress unknowns [1]. Hence,the 


following modification of Eq. 10 is required. Since 
t- 1%, %, tes? Pp (27) 


in the present case, a Mohr circle transformation must be introduced: 


7-~Mre-eMPs (28) 
and Sq. 10 is replaced by: 
= 
G = tJ (NMP) Lda (29) 


where 
! N 
cos 6 sin 8 -sin28 
™“ = sin°8 cas’ sin ad (30) 
Asigdd  - 5 sind cos 2g 


In Eq. 30, 8 is the angular coordinate of the current integration 


station on 3A. Thus, the price paid for convenience in choosing 
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P(r,6) is the slight amount of additional computation required 
to form M at each integration station, as the element boundary 


dA is swept in Eq. 29. 


2.6 Transformation-Assembly Procedure 


A complete set of stiffnesses for the near-field region can 
be obtained from a single element computation by means of a 
transformation-assembly sequence. Element number 1] in Figs. 8 
and 9 is always assumed to be oriented with edge 61 along the 
positive x-axis. Eqs. 11 and 29 are integrated with Gauss- 
Lagrange quadrature (GLQ) formulas [8} to compute Kye which is 
then assembled into a near-field substructure in accordance with 
the global node numbering scheme in Fig. 9. 

It is apparent that the stiffnesses for elements 3, 5 and 7 
can be obtained by rotation transformation of ky- For example, 
element 3 is oriented at +90 degrees from element 1. Hence, the 


transformation: 


16g tad.7 8 ee oa er (32) 


relates the locally oriented displacements 9; to the globally 


oriented displacements Gor where: 


cos9 sin®O 

-sin® cos 
cosO8 sin 
-sin8 ws 


(32) 


cos8 sin8 


(Super- diagona I) -sin@ cos (12012) 


An example for one pair of nodal displacements is illustrated in 
Fig. 9. Since strain energy is independent of any particular 
coordinate system, a global stiffness matrix for element 3 may 


be computed as follows: 


’ £ ey Fo 3 

| clea dae t. kK, H° 2 t., k, t a 
~ ~ ~ ~ 
J Introducing Eq. 31 for dp leads to 

3 ae 

YR k = ¢’ k 
HES po abst, ee 

4 : for arbitrary values of Ic° Therefore: 

7 7 

: k,° R k R 5 b= Hz (35) 
j and in a sinilar manner, 

Kook, = R kK 5 Of mM 3x/2 (36) 


The sequence of rotation transformations given by Eqs. 35 and 
36 assembly according to Fig. 9 is followed. 

A stiffness matrix for element 2 can be obtained from Ky 
by the reflection transformation iliustrated in Fig. 10. If 
element 2 is considered relative to a reflected xy axis system, 
then its stiffness is identical to k with the displacements 
oriented as indicated by the broken arrows. However, stiffness 
k, is wanted again with respect to Ic: It is evident from Fig. 


~ 


10 that for this case: 


12h. ie 


where 


5 “2 r 
A = ae rd: (38) 


{Die a . 
ym ! i (i242) 


Applying the strain energy argument again then gives immediately: 


—-— ow 


k,-ak a (39) 


Atter a has betn assembled, the procedure is completed with 
relation trarsforaations for: 

ky, & k= RKR r 

war er yg SOLAN > oF 42,6, 3m/2 (40) 
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where R is given by Eq. 32. The transformation-assembly procedure 
has been programmed internally, and the final result may be 
thought of as a single element with 16 outer boundary nodes, 24 
inner boundary nodes, and 64 total degrees of freedom. Figure 11 
illustrates the numbering convention for the assembly, which will 


be referred to as HOLEL in subsequent discussion. 


2.7 Performance Tests 

Several performance tests were conducted, in addition to 
routine checking of inversion accuracy in the H Matrix, in order 
to assess the capabilities of HOLEL. Although its primary 
function is to serve as a link between a few inner rings of 
quadrilaterals and « coarse far-field mesh, there is some interest 
in determining the solution accuracy for structures in which the 
inner boundary of HOLEL is stress-free. Errors are to be expected 
for this situation, since Che assumed surface tractions NM P & 
generally Go not vanish when P(r, 4) is computed for points Lying 
on the element boundary “A. Errors of thas typa usually appear 


as overshoots in the computed stresses, given by {i]: 


ted 


1 
ica) = Piya H &q° Bine ¢, (41) 


where r,9 represent any point in the element domain A or on the 


boundary 3A. Matrices 8B, (r,%) were computed during the formation 


i 
of ky for several points along a ray ai 4 = 27.5 deaqrees in 
element 1 (Fig. 9} for the purpose of the test. 

Figs. 3, 14 and 15 illustrate the results Sf a test series 


in which the dimensions of HOLEL are {see Fic. 12): 
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Hole Diameter Do 
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Edge Dimension W 
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The value W/D, = 4 was chosen on the dasis of a classical 
elasticity solution which indicates that far-field conditions 
are achieved -t r/R, = 4 for an open hole in an infinite plate 
(9]. The finite-element results computed from Eq. 41 are compared 
with the classical selution at 9 = 22.5 degrees [7]. 

The effect of increasing the number of GLO integration 
Stations is studied in Fig. 13, which plots Tope This stress 
component is the least sensitive to overshoot error because it 


dees not enter into the stress-free candition on the inner circular 


portion of the boundary. The results show that 3 GLQ stations 

are too few, while acceptable accuracy is obtained with 4 or more 
stations. The programming of HOLEL was subsequently fixed at 3 
5 GLQ stations basad on these results. Tho meaning of "5 stations" 
iS a&ctually: . 4 


re 


5 stations x # edges to compute S an AA = 30 


dx $ Stations to compute H in aA = 25 
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Total GLO inteqratien stations «= 45 


This amount of computing is comparable to the amount required for 
the PCRKS9 element. 
Fiqures 14 and 15 illustrate the behavior of a and Paes 


respectively. Finate-element results are shewn anly for 3 GLO 
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Stations. These results demonstrate surprisinaty that HOLEL is 


able to achieve the stress-free condition. However, some inaccuracy 


I 
| 


We 


can be observed at r/R, 1.7. The important concivsion to be 
drawn from this test series is that HOLEL can be used without 
interior elements to model a multi-fastener-hole structure for 
which, e.g., local stress distributions or K, and Ker solutions 
are sought only at one fastener hole, or at some other location. 

The second test series studied the performance of HOLEL as 
a function of the parameter W/Do: using the test problem illustrated 
in Fig. 12. Results for Te, are plotted in Fig. 16. Some performance 
degradation is observed for W/D, = 3, as the near-field stress 
gradients begin to appear in the far-field elements, which are 
incapable of following this behavior with complete accuracy. 
Considerable degradation is also seen for W/D, = 8, a result which 
might be attributed to either the absence of far-field terms in 


P(r,6), or errors in computing G and H caused by shape distortion, 


or both. One may conclude from these results that limitations 


must be placed on the edge distance and centerline spacing in 
HOLEL models of multi-fastener structures. 

The objective of the third test series was to study the 
accuracy of HOLEL for the case of combined interference-fit or 


bearing loads and panel tension. Two tests were conducted, in 


which additional loads were applied at the fastener hole boundary 


of the model shown in Fig. 12. In Fig. 17, resuits for Thy are 


plotted for the case of 1 ksi panel tension combined witn 1 ksi 


uniform pressure to represent an interference fit. A curve faired 


OS through the computed stresses extrapolates to -1.07 ksi at the 


fastener hole boundary, i.e., an error of about 7 percent. In 
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Fig. 18, results for Tt, are plotted for 1 ksi panel tension 


combined with a bearing pressure distribution given by: 


p (6) = 2f SinO 
. (42) 
P=8§000 4. , Of AS% 
Eq. 42 gives 
tpt pled= -1.95 kse (43) 


at @ = 22.5 degrees. In this case, a curve faired from 11 computed 
stresses extrapolates to the exact value at the fastener hole 
boundary. It is apparent from these results that HOLEL is 

accurate for both interference-fit and bearing loads. In the 
former case, the results could have been improved by computing 

ae at li locations instead of the 5 which had been used. Far- 
field element stresses were also checked for this test series. 

In the case of interference-fit, the computed values for cartesian 
stress in the far-field elements were Coy = 1 ksi and Cou. Oy = 0 
{to the roundoff level of the computer), with very little disturbance. 
In the case of bearing, it must be recognized that the 8,000-lb. 
bearing lcad adds to the 16,900-lb. panel tension load (1 ksi x l6- 
inch edge) for the far-field elements below the fastener hole. 

The computed stresses were in fact uy = 1 ksi for the upper 
elements, Oy = 1.5 ksi for the lower elements, and Ce oxy = 0, 


with load-transfer effects appearing in the far-field elements 


adjacent to HOLEL. 
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In the fourth test series, four rings of quadrilateral and 
PCRK59 elements were added to the structure model with the outer 
ring coupled to the circular boundary of HOLEL. Two test cases 


were analyzed. In the first case, the ring consisted entirely 


of quadrilaterals, and a cosine bearing pressure distribution 


(Eq. 42) was applied at the free boundary of the innermost ring. 


sigraaaees purer ccyare! 


Figure 19 plots the results for Tae computed at the centroids of 


peng 


one “slice" of four quadrilaterals, along a ray at 6 = 97.5 degrees. 
The faired-curve extrapolation of these four data points agrees 


well with 


2P 23 
Go - p (97.5°) =" Rg, sin (97.5°) 2 - 25,2 ksi (44) 
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In the second case, one or two groups of four quadrilaterals were 
replaced hy PCRK59 elements to simulate a fastener hole with a 
crack at 6 = 0 or two equal-length cracks at 6 = 0, 7. Several 
cases were analyzed covering the parameter ranges: 

Hole diameter/HOLEL Edge: 0.02 <D, /w < 0.05 

HOLEL edge/HOLEL diameter:W/D = 4 


Crack length/hole radius: 0.25 < a/R, < 1.0 


Table 1 compares the computed Ky values with independent classical 
solutions ,{10, 11]. The results demonstrate that the finite- 
element model is capable of providing reasonably accurate Ky solu- 
tions for quite extreme geometries (small fastener hole ina 

large panel). 


The final test series consisted of a number of demonstration 
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runs of the complete PANEL program, in which Ky and Kir solutions 
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were computed for another range of W/D, values. These results 


are presented in Section 5, following the documentation of the 


HOLEL procedure and the PANEL program. 
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HOLEL PROCEDURE 


L 

: 
Section 3 

Due to the complexity of the model, HOLEL has been programmed 

as several interrelated subroutines. Subroutine HOLEL provides 
overall control of the stiffness generation-assembly~transformation 
procedure, calling upon subroutines QHOLEL, GMTRX, HMTRX, LMTRX, 
MNMTRX, PMTRX and TRIG to execute specific computational tasks. 
In addition, HOLEL uses some of the FEABL-2 software (6] for the 


assembly process. 


| 
3.1 Structure Model and Input/Output Conventions 
A square near-field region with an inner circular boundary 
of radius Ro is modelled. ‘The numbering conventions have been 
indicated in Fig, 11. Subroutine HOLEL is invoked by: 
CALL HOLEL (COORD, THK,S,RI,RSS,ISS,B) 
where the arguments are dimensioned and defined as follows: 
COORD(12) - Vector of cartesian coordinates of the 
element corners in order Xp yr ZyeXgreees 
Za: (In the present version, the element 
is assumed to lie in the XY plane, and the 
2 coordinates need not be given any values.) 


THK - Scalar value of element thickness. 


$ (3,3) - Array vf elastic compliance constants for 


homogeneous isotropic material, i.e.: 


Cd 


RI ~ Scalar value of inner boundary radius Ro° 

RSS (2097) ,- Floating-point and fiend pein Huns of a 

ISS (2097) FEABL-2 DATA vector. These arguments must 
be equivalenced, as well as dimensioned in 


the user's MAIN program. 


B(6,3,13) - A collection of B matrices for stress analysis 


(see Subsection 3.2). 
HOLEL returns the assembled subregion stiffnesses as a variable- 
bandwidth~-stored matrix in the DATA vector RSS,ISS. The results 
may be read into another storage area with the algorithm given 
in Appendix A, or may be assembled directly into another DATA 
vector with FEABL-2 subroutine ASMSUB (see Subsection 2.2). 


Procedure HOLEL requires no input data cards. 


3.2. Reguired Subprograms and Other Features 
Procedure HOLEL requires the following additional software 


for execution: 
1. ASRL FEABL-2 subroutines ASMLTV, ORK and SETUP. 
2. IBM Scientific Subroutine Package subroutines MFSD 
and SINV. 


No external disk or tape files are required. 


‘ 1 -~P oO 
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3.3 Model Generation and Program Flow 


Figure 20 summarizes the program flow in Procedure HAOLEL. 

At entry, whatever control parameters are present in the FEABL-2 
labelled COMMON areas SIZE, BEGIN and END are saved in temporary 
storage. Subroutine SETUP is called to establish address control 
for the RSS,ISS vector for a structure model with 64 total degrees 
of freedom, consisting of the eight 12-degree-of-freedom elements 

in Fig. 9. Element interconnections are then geverated internally, 
in accordance with the numbering conventions in Fig. 8 and Fig. 1l, 
and subroutine ORK is called to compute the bund margin and complete 
address control for the assembled stiffness matrix. 

Subroutine QHOLEL is now called to compute k,, the stiffness 
matrix for element 1. This subroutine calis in turn subroutines 
HMTRX and GMTRX to compute at and G, forms the stiffnesses cTH7 tg, 
Subroutines HMTRX and GMTRX loop ove: the GLQ integration stations 
in the element domain A and on dA r-spectively, calling upon sub- 
routines TRIG, PMTRX, LMTRX and MNMTR: to compute P(r,0), L(A) 
and MN? at each station, and accumulating the product sums 
(including GLQ weighting factors) fer P'SP and (NMP) "L, Subroutine 
HNTRX calls SINV (which calls MFSD) to invert 4H. 

After ky has beer comp ‘ed, subroutine QHOLEL forms a 6x3x13 
matrix for stress analysis, ccnsisting of six matrices B(r,8) as 


given by Eq. 41, and stored in accordance with the following 


J (46) 


conventions: 


lo 


‘ 8, 
B(1,J,K) = (8; (n,4) 
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where subscript I identifies the specific B matrix, and where 


6, = X/8 5 i=4,2,--, 6 
Hh = Ro (47) 
Te= Tiyt 0.2 (distance from Ry to outer 
boundary along 6;) ; 122,3,--,6 
Thus, the B matrices refer to six locations which divide the 
ray 6 = 1/8 into five equal segments, with By located on the 
inner circular boundary and Be located on the near/far-field 
boundary. Polar stresses for these locations can be computed 


from: 


42 
Lt, Tap tre $;* STRESS (7), = Be Bayg,k)"Qtk) (48) 


where 


Q(kK) = Lo (49) 


An extra column (K=13) is appended to each B matrix for storage 
of the polar coordinates of its location. 

At this point, subroutine YOLEL regains control and executes 
the transformation-assembly process outlined in Subsection 2.6. 
In the progran, Ky is formed first, and the assembly sequence is 


programmed as: 
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4 Assembly of ky and k 
4 Rotation through "/, to k3 ana ky 
. i Assembly of k3 and k, 
Rotation through 1 to Ke and ke 
4 Assembly of ky and kg 
a Finally, the HOLEL address control parameters are transferred 
’ to the FEABL-2 labelled COMMON areas SIZESS and BEGSS, and the 
j user's control parameters are returned to areas SIZE, BEGIN and 
7 END. Control is now returned to the user's calling program, with 
HOLEL ready for assembly to the user's structure mcdel by means 
q of FRABL-2 subroutine ASMSUB. 
3.4 Procedure Status 
All of the HOLEL subroutines, together with copies of 18M 
subroutines MFSD and SINV, are maintained as a unit in the form 
j 4 of indivadually sequenced, 929-punched FORTRAN-IV source decks. 
, | A listing of Procedure HOLEL appears in Appendix B. 


Section 4 


PANEL PROGRAM 


The PANEL program is an executive program which controls 
the assembly of several substructures into a model of a skin 
tension pane] with a single open fastener hole. The PANEL program 
uses procedure HOLEL anc several] additional special-purpose proce- 
dures to generate the required substructures. Details of these 


additional procedures are discussed in Subsections 4.2 and 4.3. 


4.1 Structure Model and Input Conventions 
Figure 21 illustrates the structure which the PANEL program 
models: a rectangular panel loaded by uniform tension on its 
horizontal edges. The panel contains a single fastener hole 
eentered vertically. The fastener hole may be either centered 
er offset horisontally. The vertical edges of the panel may be 
stiffened symmetrically with intearal stiffeners, if desired. 
One or two cracks way be placed to emanate radially from the 
fastener hole. If two cracks are specified, they will be located 
180 deqrees apart, aml they may be of equal or unequal length. 
The angular position, 0, to the first crack may be varied in 
US-deyres increments from 0 to 345 degrees if there is only one 
erack, or from 0 to 145 deqrees if there are two cracks. A 


second version af the PANEL program is available for a structure 


in which onty the left cdge may be stiffened. The second version 


is otherwise similar to the first version, and will not be 
discussed separately. 

Four input data cards are re,uired for the PANEL program. 
The input data describe the specific panel ceometry and define 
several parameters which control the type of solution executed 
and the amount of information printed. Figure 22 iliustrates 
the correct format for the input cards: 


l. Panel parameters 
WIDTH 


i) 


Total width of the panel, ee 
LENGTH = Total length of the pans). 
THK = Panel thickness. 
STFFCT = Stiffener factor. 
PRESS = Tension loading, Ty" 

RI = Radius of the fastener hole. 

IOFFST = Offset indicator. 
2. Crack parameters 

A()) = Length of first crack. 

A(2) ® Length of second crack. 

IPOS (1) © Crack initial position number 

IPOS (2) = Crack final position number. 
3. Material properties 

E = Young's modulus. 

v = Poisson's ratio. 
4. Print control parameters 

RT1 = Control for optional FEABL-2 output. 


KT2 = Control for optional FARFLD output. 


KT3 = Control far optional RING output. 
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The stiffener factor is defined in terms of the panel 


dimensions (Fig. 21): 
STFFCT © wt, /Wytp (59: 


The value STSFCT = G corresponds to an unstiffened panel. Negative 
values may be used tc analyze panels with edges thinner than the 
primary structure. The offset indicator is used to control the 
type of solution desired. The following options are available; 
IOFFST = 5 - Centered fastener hole only. 
IOFFST = 1 - Fastener hole moves right. 


YOFFST = -1 - Pastener hole moves left. 


For the latter two options, solutigns are executec automatically 
beginning with the hole on canter and ending with the hole as 
close to the edge of the panel as permitted by tie finite-element 
model. If STFFCT # 0, motion of the hole is further restricted 
to avoid overlapping with a stiffener. 

A value of G should be assigned to crack length A(2) if a 
structure with only one crack is to be analyzed. Crack sizes up 
to 1.27 (RE) are permitted. The crack position numbers control 
the angular location of the first crack as follows: 

Orsni tied = (rpos(i)-4) 40 
final = {tPost2)-1) b6 han) 
602 M/i2 red. « 15 des. 


Permissible limits for the position numbers are: 


40 


{€ TPOS(A) S$ Mmay —- TPOS(1) € IPOS(Z) S Nuay (52) 


where 


Nwa, = 24 if Alz)=o 


: 53 
Nmay = 12 if Alz)#O 3) 


Solutions sre executed automatically beginning with the first 
crack in position 1 and ending with the first crack in position 
2. A Singie solution is executed if IPOS(2) = IPOS(1). A case 
Matrix is executed if ICFFST # 0 and 1P0S(2) > IPOS(1). 

Most of the routine information printed by the PANEL program 
is of no interest when production runs are executed. This informa- 
tion may be deleted from the output by assigning each of the 
contro] parame cers KTI1, KT2, KT3 values equal te the FORTRAN 
unit nuwaber for the line printer at the user's computing facility.* 
Any Other value permits full outout by the associated software. 
Full output is recommended for initial testing of the pvroeqram at 
anew facility. Output from FARPLD and RING shovld be allowed 


whenever a new range of panel dimensions is tried. 


4.2 Reaeired fubpresyr ars ava Other Features 
et ee, epee MP es ae ere) nee Mle See Re a cee ners. SA emt) Ne ee eee Set a 


The PANEL program requires the following additional software 


for execution: 


me ie a ee ee ee 


*The same value used in a print instruction, e.q., WRITE (6, 1000} 
a,B,C. The FORTRAN unit number in this case is 6. 


Reet AS ren ei mer eee gine 


l. ASRL FEABL-2 subroutines ASMLTV, ASMSUB, BCON, FACT, 
ORK, QBACK, SETUP, SIMULQ, STACON and XTRACT. 
2. IBM Scientific Subroutine Package subroutines MFSD 
and SINV., 
3. ASRL procec'res FARFLD, HOLEL and RING, with included 
Subroutines. 
4. ASRL elements PCRK59 and QUAD4. 
No external disk or tape files are required. The program must be 
able to communicate with the user’s facility card reader und line 
printer, by means of two instructions near the heginning of the 
program: 
KR = Card reader FORTRAN unit number. 
KW = Line printer FORTRAN unit number. 


The progran is supplied with IBM-standard values KR=5 and KwWe6. 


4.3 Model Generation and Program Flow 

The panel structure model is created via several levels of 
substructuring to minimize repeated processing of unwanted degrees 
of freedom. Figure 23 illustrates the general hierarchy of the 
finite-element model. Procedures PARFLD and RING generate the 
two major components which are finally assembled to form the 
conplete structure, 

Procedure FARFPLD creates the far-field region, consisting 
of upper and lower rectangular substructures generated by 
procedure LUG,* and a center-zsone substructure generated by 


We Nae de a ee Ne Lerma nee enema 


*pifferent from the attachment luq procedure reported in Ref. i. 
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procedure CZONE. Procedure LUG assembles a regular mesh of 

QUAD4 elements and eliminates all except the upper-edge and 
lower-edge nodes. Procedure CZONE assembles HOLEL together with 
right and left portions modelled with QUAD4 elements, and eliminates 


all but the upper- and lower-edge nodes and the inner circular 


boundary uodes. Procedure FARFLD assembles the LUG and CZONE 
components and executes additional Gauss elimination to produce 
one of two major substructures; 

1. "Panel-and-Hole", in which only the nodes along the 
top and bottom edges of the panel and the nodes along 
the inner circular boundary remain as boundary nodes 
in the statically condensed structure. Displacement 
restraints are applied at the bottom edge and nodal 
forces are applied at the top edge to represent 
uniform tension. 

2. "The Chcshire Cat": a panel-and-hole with top and 
bottom edges eliminated by static condensation (nothing 
remains but the smile). 

Only the Cheshire Cat option is used by the PANEL program. Number- 
ing conventions for the FARFLD compcnents are illustrated in Figs. 
24, 25 and 26. 

Procedure RING creates an inner cracked riny for assembly 
inside the Cheshire Cat. The ring consists of two 150-degree 
arcs and two 30-degree segments which contain QUAD4 and PCRK59 


elements. The arcs are assembled by procedure ARC4, which also 


per ee Gt aOR ys 7 
Rey Ae Ne bday Medan stash cua ddate « huecc seabed mee thee A cdettnbe 
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removes all interior nodes by static condensation. Procedure 
RING assigns the angular locations of the components to obtain 
the required crack positions, and determines where each PCRK59 
element is to be placed in the 30-degree segments according to 
the crack sizes A({1), A(2). Finally, procedure RING assembies 
and statically condenses the cracked ring to produce one of two 
Major substructures: 
1. Ring with all but inner and outer circular boundary 
nodes eliminated. 
2. Ring with all but outer circular boundary nodes 
eliminated. 
Retention of the inner boundary nodes is useful for cases in which 
bearing or interference-fit loads are to be applied at the fastener 
hole. The outer boundary nodes must be retained for coupling 
with the Cheshire Cat. The PANEL program uses only the second 
option. Numbering conventions for the RING components are 
illustrated in Figs. 27 through 30. Procedure RING always 
generates a ring with an cuter/inner diameter ratio of 2.52 to 
maintain shape conformity for the QUAD4 elements. 

Figure 31 illustrates the executive flow in the PANEL program. 
After the problem input data have been read and printed and some 
Parameters have been calculated for control of the various 
procedures, two major loops appear. The outer loop begins with 
the creation of a Cheshire Cat, a step which must be repeated 
each time the fastener hole is offset to a new position. This 


is followed by address control and interconnection generation for 
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the final structure, which will be azsembled from the Cheshire Cat 
{(super-element 1) and the cracked ring (super-element 2). The 
interconnection algorithm is actually completed inside the inner 
loop to zllow the band-margin computations (Subroutine ORK) to 
erase the K and q data from the previous pass. The remainder of 
the inner loop consists of assembly of the Cheshire Cat, generation 
and assembly of the cracked ring, global solutien (subroutines 

FACT and SIMULQ}, a back-substitution to obtain qr, for the cracked 
ring substructure and finally, extraction of the PCRK59 displace- 


and K, The ends of the 


I I° 
inner and outer loops are governed by incrementation of the crack 


ments from qr, and computation of K 


angle and hole offset, respectively, and logical checks to determine 


whether the prescribed ranges of these parameters have been swept. 


4.4 Output Conventions and Error Messages 


Figure 32 illustrates a sample output from the PANEI. program 
(version 2, lett side stiffened). The output from Version 1 is 
Similar, The heading identifies the program version and repeats 
the user's input data. Below the material properties data appears 


a table of Ky and K solutions tor one or two cracks, together 


II 
with their angular positions. The Ky and K values are NASA/ASTM 


II 
standard stress intensities in unics of psi Yin., assuming that 
the input data was specified in corresponding engineering units 
of psi for loading and Young's modulus and inches for dimensions. 
The sample output is an example of the production information 


obtained by exercising the three options for deletion of debugging 


output. 
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Besides the FEABL-2 software error messages [6], the follow- 
ing diagrostics may result from the PANEL programs. In procedure 
RING, the lenath specified for each crack is checked to insure 
that the crack tip does not extend beyond the average radius of 
the outermost ring of quadrilaterais, i.e.: 

3 
Ry + AIT) S$ 5 CRG H+ 2.52) 5 THh,2 (54) 

where R, is the inner radius cf the outermost ring. If Eq. 54 
is violated, a message is printed and program execution is terminated. 
T.is condition is somewhat conservative, since it restricts the 


crack tip to a position mid-way in the PCRK59 element. Eq. 54 


may be replaced by: 
Rpt AlT)S d3Re + 0.71 252Rz (55) 


to permit the crack tip to approach somewhat closer to the outer 
bour.Jary. In procedure CZONE, a check is made to insure that the 
fastener hole offset is within allowable limits. This is done by 
monitoring the node number of the fictitious center reference node 
(see Fig. 25). Excessive offset causes an error message and 


program termination. 


4.5 Program Status 


Both versions of the PANEL program have been exercised 


successfully for all analysis options, and are maintained as 
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sequenced, 029-punched FORTRAN-IV source decks. Either version 


requires 260 to 300 KBYTES (65,000 to 75,000,, words, or 


10 10 


176,750, to 222,370, words) of core memory and approximately 


8 
0.8 to 1.0 CPU minute, depending upon the ranges of 

the hole offset and crack angle parameters. Storage and time 
Statistics are based on runs made on an IBM S-370/168 machine, 
using the IBM FORTRAN-Gi and FORTRAN-H compilers. Version 1 


(symnetrical edge stiffeners) is listed in Appendix C. Version 


2 (left edge stiffener only) is listed in Appendix D. 
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Section 5 


DEMONSTRATION EXAMPLES 


A number of example analyses have been run to verify the 
PANEL program code and simultaneously to explore the accuracy of 
the analysis for a range of panel and crack dimensions wider 
than considered in the HOLEL performance tests. The solutions 
given here focus mainly on panels with centered fastener holes 
having cracks at @=0,1 because there exist no independent solu- 
tions with which to compare other configurations. 

Figure 33 illustrates the coarsest model possible to generate 
from the PANEL program: lugs consisting of two elements each and 
a center-zone composed only of HOLEL. Also, the dimensions 
chosen are such that the fastener hole is no longer small compared 
to the panel, and since the outer/inner diameter ratio of the 


cracked ring is 2.52, the HOLEL shape parameter in this case is: 
W/D, = 4/2.52 = 1.59 


a point well outside the range studied in the HOLEL performance 


tests. The butterfly plots for Ky and K shown in Fig. 33 


It 
illustrate an error effect caused by proximity of the displacement 
boundary conditions to the region of interest. In the present 
case, the restraints are only two elements away from the “action" 


{a lug QUAD4 and the HOLEL) when either crack lies below the 
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horizontal. The error is most pronounced for Kyy at 8 = 225, 315 
degrees. The error would be more pronounced as 6 + 270 degrees, 

except that the solution tends rapidly to zero in this region. 

A similar but less-pronounced effect can be observed in the plot 


for K,. Comparison of K. for 6 = 0, 180 degrees with an indepen- 


I I 
dent solution [10, 11] provides a more welcome result. Indeed, 
it is quite surprising that this very crude finite-element model 
is able to faithfully reproduce the classical solution. The 
general conclusion to be drawn from this test is that only the 
upper half of the butterfly plot may be trusted when running 
models with very few elements in the FARFLD substructure. 

Figures 34 and 35 illustrate two runs in which the panel 
and crack dimensions have been varied to explore the effect of 
W/D,° Also, the aspect ratios of the elements in the LUG sub- 
Structure were changed to allow more elements between the restrained 
bottom edge and the center-zone. The latter modification has 
eliminated the restraint error effect. Table 2 summarizes the 
comparisons of results from Figs. 33 through 35 with the independent 


solution, showing that reasonable accuracy has been achieved over: 
1.59 $ W/D, @ 6.35 


The greater error for the middle case remains an unexplained 
anomaly. 


Figure 36 illustrates the K, and K solutions obtained for 


T If 
a panel of the same dimensions, hole radius, ete. scown in Fiq. 


35, but with the hole offset to the maximum amount permitted in 
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the model, a distance of 1.5 inches. Increases of the order 
of 5 to 10 percent are observed for Ki, with the greatest 
increase at crack tips which are located nearest to the edge 
of the panel. 

Figures 37 through 39 present results for some additional 
studies of the basic panel shown in Fig. 35 to illustrate the 
stiffening options. A stiffener factor of 0.5 was used for 
these analyses, i.e., the added cross section area of each 
stiffener was 50 percent of the panel cross section area. 
Figure 37 plots K 


and Ky solutions for a symmetrically 


I I 
stiffened panel. Results are shown for the fastener hole 

centered and offset by 1 inch. There is very little difference 
between the two solutions. Solutions for a panel with left-edge 
stiffener and a centered fastener hole are plotted in Fig. 38. 

The results for Ky and Kry appear to be symmetrical. In fact, 

very little difference can be observed between the unstiffened, 
symmetrically stiffened and asymmetrically stiffened panels (compare 
Figs. 35, 37 and 38). Some diffe~ence is noted when the fastener 
hole is offset 1.5 inches to the right, away from the stiffened 
edge. Figure 39 presents the results for this case, which are 
almost identical to the results for an unstiffened panel (compare 
with Fig. 36). The tentative conclusion from these results is 

that edge stiffeners and moderate offsets do not appreciably affect 
K, and Kyy for cracks at the fastener hole, at least for uniform 


T 
tension loading and when the hole and crack are not extremely close 
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to an edge or stiffener. 

Figure 40 summarizes the dimensional information for a panel 
model used to study a broad range of a/R, ratios with W/D, fixed 
at 1.59. Table 3 summarizes the results for a series of cases 
with equai-length cracks at 6=0,n. These may be compared directly 
with handbook data because the dimensions of the model correspond 
precisely with a published curve. The results are seen to be 


quite reasonable over the entire range: 
0.05 ¢ a/R, € 1.25 


Indeed some of the error indicated in the table may be attributed 
to inaccuracy in reading the handbook chart. A number of other 
cases for a single crack at 6=0 and for two unequal-length cracks 
are compared with the foregoing results in Table 4. Complete 
butterfly plots of the Ky and Kry data computed from these runs 


will appear in a later report. 
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Section 6 


CONCLUSIONS 


This report has traced the development of a parametric 
finite-element analysis program for computation of Mode I and 
Mode II stress intensity factors in stiffened and unstiffened 
panels with cracks emanating from centered and offset fastener 
holes. Early attempts to model the structure entirely with 
conventional elements, except for the crack-tip element and a 
few mesh-expander elements, proved to be unsuccessful. Test 
runs of these types of finite-element models quickly demonstrated 
that asymmetric grading of the mesh was forced by the need to 
accommodate an offset hole and at the same time to keep the size 
of the mode} within economic bounds. The asymmetric mesh was 
found to give extremely poor results for compute ! displacements 
in regions of low stress gradient, and was therefore abandoned. 

The assumed-stress hybrid method was called upon once again, 
this time to provide a special element which could handle the 
transition from circular geometry near the fastener hole to the 
cartesian geometry natural to the rest of the panel. As finally 
developed, the hybrid element occupied a 45-degree sector between 
its inner (circular) and outer boundary. Airy stress functions 
from classical elasticity solutions near a hole in an infinite 


plate were used to provide the assumed stress field. The new 


element allowed a much cleaner, more economical finite~element 
mesh to be designed, while at the same time providing accurate 
computations near the fastener hole. Performance tests have 


shown that the element is capable of occupying a stress-free 


fastener hole boundary directly, and of accepting interference- 
fit and cosine bearing loading on the fastener hole surface. 
Additional tests demonstrated that conventional quadrilateral 
elements could be ceupled in rings inside the inner boundary of 
the new hybrid element, with no loss of accuracy. Some tests were 
also conducted with hybrid crack-cantaining elements replacing 2 
some »f the quadrilaterals. These tests demonstrated that stress 
intensity factors could be computed to within a few percent of 

the values given by well-established independent selutions based 

on classical methods. In the final phase of the project, 4 

parametric program was developed and verified for analysis of 

tension panels if the various configurations mentioned above. 

During the verification tests, the hybrid fastener hole-ring 

combination was subjected te a variety of dimensional and shape 

parameters to further extend the range of measured performance. 

The results of these tests have demonstrated that the panel 

pragram is cabable of computing stress intensity factars to 

within § percent or better for fastener hole ami crack sizes found 

in current Nieecawes: 


The many example results presented im this report are still 


A rather limited data base, when compared with the number of atress 


intensity factor solutions needed for a comprehensive designer's 
handbook. Some additional dats, which were generated in the 
final verification tests but not included in this report, will 
appear in a later report in this series. However, further verifica- 
tion tests are still required to broaden the range of applicability 
ef the panel program. One concern which has not yet been answered 
is how close a row of fastener holes may be spaced, or how near 
may & fastener hole approach the edge of a panel, before the 
finite-element model experiences unacceptable degradation of 
accuracy. The data generated thus far seem to indicate that there 
with be a limit in this respect, and that the limit may be severe. 
Addition of mid-edge nodes to the edges which span between the 
hybrid fastener hole element's inner and outer boundaries may 
berwit closer spacing, but will also require additional terms 
hh the assupad stress field. 

| Another continuing concern is associated with the verification 
process itself. The combination of crack-containiag and fastener 


hole elements gives the capability to compute K, and K for so 


i It 

many varied configurations that the numerical analyst finds hia- 

seif sailing in uncharted waters. At the present time, parametric 
codes like the panel program cai only be calibrated against classical 
solutiona for Ky at ane or two data points, while they may compute 


as many as 100 data points each for K, aud Kegs 
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TABLE 1] 


COMPARISON OF CLASSICAL AND FINITE~ELEMENT 
FRACTURE MECHANICS SOLUTIONS 


Number of Crack Size, Hole Radius, K. (psi/in.) %, 
Cracks a (in.) Ry (in.) a/R, Bact FEM Error 
2 0.04 0.08 0.5 650 630 3.0 
1 0.04 0.10 0.4 660 675 2.4 
2 0.08 0.10 0.8 792 771 2.7 
2 0.05 0.20 0.25 910 910 oe 
2 0.20 0.20 1.0 1150 1170 i.9 
TABLE 2 
COMPARISON OP CLASSICAL AND FINITE-ELEMENT 
RESULTS FROM THREE DEMONSTRATION EXAMPLSS 
(Fastener hole with 2 cracks} 
Crack Hole HOLEL X, (psiein) % 
Case Fiq. Size {in.) Rad. (in.) a/R, W/D,, ERact FEM Error 
i 33 q.3 0.5 0.6 2.59 12.9  exao 1.0 
2 34 0.125 0.25 O.s 3.47 pai ii)t4 4.1 
J 35 0.125 O.285 1.0 &.35 B96 §Ha a 7 
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TABLE 3 


PERFORMANCE OF PANEL PROGRAM AS A 
FUNCTION OF a/R, RATIO 


(Structure model in Fig. 40; cracks at @=0,17.) 


a(in.) and K_ (psiVin.) g 
a/R, Exdct* FEM Error 
0.05 1144 1169 1.4 
0.1 1673 1613 3.6 
0.2 2037 2052 0.7 
0.4 2474 2445 1.2 
0.6 2735 2715 0.7 
0.8 2973 2897 2.6 
1.0 3257 3099 4.9 
1.2 3521 3491 3.4 
1.25 3588 3500 2.5 


*Ref. 11, p. 19.4, curve marked h/b=2, R/b=0.25. 


TABLE 4 


ADDITIONAL TEST OF a/R, RATIO 


(Structure model in Fig. 40; one crack or two unequal cracks. 


a(in-) and a/R, K. (psiYin.) at =0 % 
8=0 9=1 5 FEM Exact (from Table 3) DXifference 
0.05 ~-- 1155 1144 1.0 
0.1 --- 1600 1673 4.4 
0.2 --- 2000 2037 1.8 
0.4 --- 2302 2474 6.9 
0.6 --- 2466 2735 9.8 
0.8 wo 2522 2973 15.2 
1.0 --- 2602 3257 20.) 
1.2 --- 2756 3521 7 
1.25 --- 2813 588 21.6 
0.05 1.25 1447 1144 26.5 
1.25 0.05 2824 3588 21.3 
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Appendix A 


ALGORITHM FOR TRANSFERRING A SUBSTRUCTURE 


TO RESERVE STORAGE 


Assume that a standari TLABL-2 assembly has been executed, 
and that the assembied structure has been statically condensed 


with subroutine STACON. The objective is two transfer the con- 
Ge eat eels AAC) = 
ROB and force Qn from the FEABL DATA 


vector (desicnated by RNANE,INAME) to reserve storage. Assume 


densed stiffness matrix K 


: she? x — Gs 2 fo). 
further that x! is almost fully populated. Therefore, K “Fo is 
~EB oe wet xe - BB 


to be inflated to a fully populated, ilower-triangle-vector-stored 
matrix during the transfer gerocess. For this pursese, vectors 
STORK and STORQ are dimensioned in the MAIN program ta at least 
NiNel}s2 and N words respectively. where N is the total number 
ef boundary €uncondensed) degrees of Freedom remaining. Vector 


aa pane ye 
RA >» white STOR receives Oo: . 
“5 


a: one 


STORE will recerave K 
Ocher var.ables appearic,; ino the ataeriphm play the fallow he 
retas, ISERG as the Farab uncendensed dearee af Freedem, 4.0., 


L4ERO * NIO¢1. To oand J are leep indices wach cantrol progress 


5 fo) te) fs : 
through K meee: ASS an the RATA wecter: | {or the .2twe ant ! 
RR oa 


for the ealumns. fF « rew is enenunteread far which the band mercin 


{ct 


of Kop 


lies to the right of column TEESO, then suize loading zeros 


‘must be inserted in the corresponding raw in STORK. LNZCOL is 


used to compute and compare the band margin. Ti is used ta bald 


od 


the variable-bandwidth address infurmation required to lncate 
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stiffnesses in the DATA vector belonging to ROW I. M, N and MN 
are used to trace the lower-triangle address in vector STORK. 
IKOUNT, ILNZ and IQ are FEABL address control parameters located 
in the BEGIN labelled COMMON area. NDT is the total number of 
degrees of freedom in the assembled structure (boundary plus 
interior), available in he SIZE labelled COMMON area. 

The transfer algorithm begins immediately after subroutine 


STACON has been executed: 


IZERO = NID+1 
M=0 
DO 20 I=IZERO, NDT 
II=. NAME (IKOUNT+I-1) 
M=M+1 
STORQ (M) =RNAME (TQ+I-1) 
N=0 
DO 10 J=I4ZERO,I 
N=N+1 
MN=(M* (M=1) ) /24+N 
STORK (MN) =0.0 
LNZCOL=INAME (ILNZ+I-1) 
IF (LNZCOL-J) 5,5,10 

5  STORK(MN) = RWAME(II+J) 

19 CONTINUE 


20 CONTINUE 
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